GitHub

Example: Reachability problem solved by Lazy ellipsoid abstraction.

using StaticArrays, Plots
using JuMP, Clarabel

import Random
Random.seed!(0)

using Dionysos
const DI = Dionysos
const UT = DI.Utils
const ST = DI.System
const PR = DI.Problem
const OP = DI.Optim
const AB = OP.Abstraction

using Symbolics

include(joinpath(dirname(dirname(pathof(Dionysos))), "problems", "non_linear.jl"))
Main.NonLinear

First example

concrete_problem = NonLinear.problem()
concrete_system = concrete_problem.system
Dionysos.System.SymbolicSystem(Main.NonLinear.var"#4#5"(), Symbolics.Num[1.1px - 0.2py + vx - 5.0e-5(py^3), 0.2px + 1.1py + vy + 5.0e-5(px^3)], 1.0, 2, 2, 2, Symbolics.Num[px, py], Symbolics.Num[vx, vy], Symbolics.Num[wx, wy], IntervalArithmetic.Interval{Float64}[[-1.0, 1.0]_com, [-1.0, 1.0]_com], IntervalArithmetic.Interval{Float64}[[-20.0, 20.0]_com, [-20.0, 20.0]_com], IntervalArithmetic.Interval{Float64}[[0.0, 0.0]_com, [0.0, 0.0]_com], Dionysos.Utils.HyperRectangle{2, Float64}([-20.0, -20.0], [20.0, 20.0]), Dionysos.Utils.HyperRectangle{2, Float64}([-10.0, -10.0], [10.0, 10.0]), Dionysos.Utils.HyperRectangle{2, Float64}([0.0, 0.0], [0.0, 0.0]), Dionysos.Utils.Ellipsoid{2, Float64, StaticArraysCore.SMatrix{2, 2, Float64, 4}, StaticArraysCore.SVector{2, Float64}}[Dionysos.Utils.Ellipsoid{2, Float64, StaticArraysCore.SMatrix{2, 2, Float64, 4}, StaticArraysCore.SVector{2, Float64}}([0.03333333333333333 0.0; 0.0 0.03333333333333333], [0.0, 0.0])], Main.NonLinear.var"#f_eval#system##0"{Float64, Float64}(1.0, 5.0e-5), Main.NonLinear.var"#f_backward_eval#system##1"{Float64, Float64}(1.0, 5.0e-5), [[0.1 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.1]], [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0])

Optimizer's parameters

sdp_opt = optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true)

maxδx = 100
maxδu = 10 * 2
λ = 0.01
k1 = 1
k2 = 1
RRTstar = false
continues = false
maxIter = 300

optimizer = MOI.instantiate(AB.LazyEllipsoidsAbstraction.Optimizer)
AB.LazyEllipsoidsAbstraction.set_optimizer!(
    optimizer,
    concrete_problem,
    sdp_opt,
    maxδx,
    maxδu,
    λ,
    k1,
    k2,
    RRTstar,
    continues,
    maxIter,
)

Build the state feedback abstraction and solve the optimal control problem using RRT algorithm.

MOI.optimize!(optimizer)
Iterations2Go:	300
	Closest Dist: 28.284271247461902
Iterations2Go:	299
	Closest Dist: 28.284271247461902
Iterations2Go:	298
	Closest Dist: 28.284271247461902
Iterations2Go:	297
	Closest Dist: 28.284271247461902
Iterations2Go:	296
	Closest Dist: 28.284271247461902
Iterations2Go:	295
	Closest Dist: 28.284271247461902
Iterations2Go:	294
	Closest Dist: 28.284271247461902
Iterations2Go:	293
	Closest Dist: 28.284271247461902
Iterations2Go:	292
	Closest Dist: 28.284271247461902
Iterations2Go:	291
	Closest Dist: 28.284271247461902
Iterations2Go:	290
	Closest Dist: 28.284271247461902
Iterations2Go:	289
	Closest Dist: 28.284271247461902
Iterations2Go:	288
	Closest Dist: 28.284271247461902
Iterations2Go:	287
	Closest Dist: 28.284271247461902
Iterations2Go:	286
	Closest Dist: 28.284271247461902
Iterations2Go:	285
	Closest Dist: 24.494026507183737
Iterations2Go:	284
	Closest Dist: 12.792689040381546
Iterations2Go:	283
	Closest Dist: 12.792689040381546
Iterations2Go:	282
	Closest Dist: 12.792689040381546
Iterations2Go:	281
	Closest Dist: 12.792689040381546
Iterations2Go:	280
	Closest Dist: 12.792689040381546
Iterations2Go:	279
	Closest Dist: 12.792689040381546
Iterations2Go:	278
	Closest Dist: 12.792689040381546
Iterations2Go:	277
	Closest Dist: 12.792689040381546
Iterations2Go:	276
	Closest Dist: 6.244347726304238
Path cost from EI : 1040.8819123608534

Get the results

abstract_system = MOI.get(optimizer, MOI.RawOptimizerAttribute("abstract_system"))
abstract_problem = MOI.get(optimizer, MOI.RawOptimizerAttribute("abstract_problem"))
abstract_controller = MOI.get(optimizer, MOI.RawOptimizerAttribute("abstract_controller"))
concrete_controller = MOI.get(optimizer, MOI.RawOptimizerAttribute("concrete_controller"))
abstract_lyap_fun = MOI.get(optimizer, MOI.RawOptimizerAttribute("abstract_lyap_fun"))
concrete_lyap_fun = MOI.get(optimizer, MOI.RawOptimizerAttribute("concrete_lyap_fun"));

Simulation

We define the cost and stopping criteria for a simulation

cost_eval(x, u) = UT.function_value(concrete_problem.transition_cost, x, u)
reached(x) = x ∈ concrete_problem.target_set
nstep = typeof(concrete_problem.time) == PR.Infinity ? 100 : concrete_problem.time; # max num of steps

We simulate the closed loop trajectory

x0 = concrete_problem.initial_set.c
x_traj, u_traj = ST.get_closed_loop_trajectory(
    concrete_system,
    concrete_controller,
    x0,
    nstep;
    stopping = reached,
    f_map_override = (x, u) -> concrete_system.f_eval(x, u, [0, 0]),
)
c_traj, cost_true = ST.get_cost_trajectory(x_traj, u_traj, cost_eval)
cost_bound = concrete_lyap_fun(x0);
println("Goal set reached")
println("Guaranteed cost:\t $(cost_bound)")
println("True cost:\t\t $(cost_true)")
Goal set reached
Guaranteed cost:	 1040.8819123608534
True cost:		 818.7712347342838

Display the results

Display the specifications and domains

fig = plot(;
    aspect_ratio = :equal,
    xtickfontsize = 10,
    ytickfontsize = 10,
    guidefontsize = 16,
    titlefontsize = 14,
    label = false,
);
xlabel!("\$x_1\$");
ylabel!("\$x_2\$");
title!("Specifictions and domains");

#Display the concrete domain
plot!(concrete_system.X; color = :grey, opacity = 0.5, label = false);

#Display the abstract domain
plot!(abstract_system; with_arrows = false, cost = false, label = false);

#Display the concrete specifications
plot!(concrete_problem.initial_set; color = :green, label = false);
plot!(concrete_problem.target_set; color = :red, label = false)
Example block output

Display the abstraction

fig = plot(;
    aspect_ratio = :equal,
    xtickfontsize = 10,
    ytickfontsize = 10,
    guidefontsize = 16,
    titlefontsize = 14,
);
title!("Abstractions");
plot!(abstract_system; with_arrows = true)
Example block output

Display the Lyapunov function and the trajectory

fig = plot(;
    aspect_ratio = :equal,
    xtickfontsize = 10,
    ytickfontsize = 10,
    guidefontsize = 16,
    titlefontsize = 14,
);
xlabel!("\$x_1\$");
ylabel!("\$x_2\$");
title!("Trajectory and Lyapunov-like Fun.");

for obs in concrete_system.obstacles
    plot!(obs; color = :black)
end
plot!(abstract_system; with_arrows = false, cost = true);
plot!(x_traj; color = :black)
Example block output

This page was generated using Literate.jl.